Sea ice thickness

In this notebook, we’ll show more detailed analysis and visualization of ICESat2 sea ice thickness measurements and the PIOMAS (Pan-Arctic Ice Ocean Modeling and Assimilation System) modelled sea ice thickness data. We’ll also highlight interactive mapping using the hvplot and holoviews package.

import xarray as xr # For working with gridded climate data 
import pandas as pd
import numpy as np
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.misc_utils import restrictRegionally # Region restriction 
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
import holoviews as hv
from textwrap import wrap
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore') 

1) Read in the data

We’ve included a function in the read_data_utils module that reads in the book dataset from the public google storage bucketas an xarray.Dataset object.

print(read_book_data.__doc__) # Print docstring
 Read in data for ICESat2 jupyter book. 
    If the file does not already exist on the user's local drive, it is downloaded from the books google storage bucket (https://console.cloud.google.com/storage/browser/is2-pso-seaice)
    The netcdf file is then read in as an xr.Dataset object 
    
    Args: 
        None
    Returns: 
        book_ds (xr.Dataset): data 
    
    
book_ds = read_book_data() # Read/download the data 
print(book_ds) # Show dataset
<xarray.Dataset>
Dimensions:                  (time: 30, y: 448, x: 304)
Coordinates:
  * time                     (time) datetime64[ns] 2018-11-01 ... 2021-04-01
    longitude                (y, x) float32 ...
    latitude                 (y, x) float32 ...
    region_mask              (y, x) uint8 ...
Dimensions without coordinates: y, x
Data variables: (12/18)
    freeboard                (time, y, x) float32 ...
    ice_density              (time, y, x) float32 ...
    ice_thickness            (time, y, x) float32 ...
    ice_thickness_smoothed   (time, y, x) float64 ...
    ice_thickness_unc        (time, y, x) float32 ...
    ice_type                 (time, y, x) float32 ...
    ...                       ...
    piomas_ice_thickness     (time, y, x) float32 ...
    t2m                      (time, y, x) float32 ...
    msdwlwrf                 (time, y, x) float32 ...
    drifts_uT                (time, y, x) float32 ...
    drifts_vT                (time, y, x) float32 ...
    drifts_magnitude         (time, y, x) float32 ...
Attributes:
    description:    Data used in ICESat-2 jupyter book
    note:           See individual data variables for references
    creation date:  2021-09-20

Restrict data to a given region

We’ve built a region mask for the Arctic into the dataset used for this book; see the data wrangling notebook for more information. The region mask is included as a coordinate in the dataset, allowing us to easily access the different regions. You can restrict the data by selecting keys corresponding to regions of interest.

We are often interested in a region called the Inner Arctic, which is defined as the combined area of the Central Arctic, Beaufort Sea, Chukchi Sea, E Siberian Sea and the Laptev Sea. Figure from Petty et al., (2020):

Inner Arctic region map

Each region is defined by a specific integer key value. We can select the region/s defined by the key/s by using the restrictRegionally function, defined in the misc_utils module. Here, we’ll restrict the data to the Inner Arctic by selected the key values [10,11,12,13,15], corresponding to the five regions described above.

regions = pd.DataFrame({"labels":book_ds.region_mask.attrs["labels"]}, index=book_ds.region_mask.attrs["keys"])
print(regions)
                                     labels
0                     Lakes, extended coast
1                         non-region oceans
2                  Sea of Okhotsk and Japan
3                                Bering Sea
4                                Hudson Bay
5                      Gulf of St. Lawrence
6   Baffin Bay, Davis Strait & Labrador Sea
7                             Greenland Sea
8                              Barents Seas
9                                  Kara Sea
10                               Laptev Sea
11                        East Siberian Sea
12                              Chukchi Sea
13                             Beaufort Sea
14                     Canadian Archipelago
15                             Arctic Ocean
20                                     Land
21                                    Coast

The regions defining the Inner Arctic are associated with the key values 10, 11, 12, 13, and 15. Here, we’ll use the restrictRegionally function defined in the misc_utils module to mask out data outside of that region.

innerArcticKeys = [10,11,12,13,15] #Inner Arctic
book_ds_innerArctic = restrictRegionally(book_ds, innerArcticKeys)
Regions selected: Inner Arctic

On the map below, you can see how the data we’ve called the function on has been cut to show only gridcells in the Inner Arctic!

date = "Feb 2020"

fig, axes = plt.subplots(1, 2, figsize=(5,5), subplot_kw={'projection':ccrs.NorthPolarStereo(central_longitude=-45)})
im1 = book_ds["ice_thickness_smoothed"].sel(time=date)[0].plot.pcolormesh(ax=axes[0], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
im2 = book_ds_innerArctic["ice_thickness_smoothed"].sel(time=date)[0].plot.pcolormesh(ax=axes[1], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
plt.colorbar(im1, cax=fig.add_axes([0.93, 0.25, 0.025, 0.5]), extend='both', label="sea ice thickness (m)")

axes[0].set_title("\n".join(wrap("Sea ice thickness before restricting, "+date, 32)), fontsize=8)
axes[1].set_title("\n".join(wrap("Inner Arctic sea ice thickness, "+date, 32)), fontsize=8)
axes[0].coastlines()
axes[1].coastlines()

plt.show()
_images/sea_ice_thickness_11_0.png

Set years over which to perform the analysis

Setting years = [2018,2019,2020] indicates that you want to perform the analyses in the notebook for three winter season: winter 2018-19, 2019-20, and 2020-21

years = [2018,2019,2020]

3) ICESat-2 interpolated sea ice thickness maps

Here, we’ll overlay sea ice thickness data over a map of the Arctic. To generate the map, we’ll use interactive plotting functions based off the hvplot python package. To see more information on the functions used, click on the plus buttons to see the function documentation, or go into the plotting_utils module to see the code.

print(interactiveArcticMaps.__doc__) # Print docstring 
 Generative one or more interactive maps 
    Using the argument "slide", the user can set whether each map should be displayed together, or displayed in the form of a slider 
    To show each map together (no slider), set slider=False
    
    Args: 
        da (xr.Dataset or xr.DataArray): data to restrict by time; must contain "time" as a coordinate 
        clabel (str, optional): colorbar label (default to "long_name" and "units" if given in attributes of da)
        cmap (str, optional): matplotlib colormap to use (default to "viridis")
        colorbar (bool, optional): show colorbar? (default to True)
        vmin (float, optional): minimum on colorbar (default to 1st percentile)
        vmax (float, optional): maximum on colorbar (default to 99th percentile)
        title (str, optional): main title to give plot (default to no title)
        ylim (tuple, optional): limits of yaxis in the form min latitude, max latitude (default to (60,90))
        frame_width (int, optional): width of frame. sets figure size of each map (default to 250)
        slider (bool, optional): if da has more than one time coordinate, display maps with a slider? (default to True)
        cols (int, optional): how many columns to show before wrapping, if da has more than one time coordinate (default to 3)
    
    Returns: 
        pl (Holoviews map)
    
    
pl = interactiveArcticMaps(book_ds["ice_thickness_smoothed"], colorbar=True, slider=True, frame_width=500)
pl
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
/var/folders/8v/cr06mz0n3bjd5jr12_9v6n200000gn/T/ipykernel_35719/2060953563.py in <module>
      1 pl = interactiveArcticMaps(book_ds["ice_thickness_smoothed"], colorbar=True, slider=True, frame_width=500)
----> 2 pl

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/displayhook.py in __call__(self, result)
    260             self.start_displayhook()
    261             self.write_output_prompt()
--> 262             format_dict, md_dict = self.compute_format_data(result)
    263             self.update_user_ns(result)
    264             self.fill_exec_result(result)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    149 
    150         """
--> 151         return self.shell.display_formatter.format(result)
    152 
    153     # This can be set to True by the write_output_prompt method in a subclass

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    148             return {}, {}
    149 
--> 150         format_dict, md_dict = self.mimebundle_formatter(obj, include=include, exclude=exclude)
    151 
    152         if format_dict or md_dict:

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in __call__(self, obj, include, exclude)
    968 
    969             if method is not None:
--> 970                 return method(include=include, exclude=exclude)
    971             return None
    972         else:

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/core/dimension.py in _repr_mimebundle_(self, include, exclude)
   1315         combined and returned.
   1316         """
-> 1317         return Store.render(self)
   1318 
   1319 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/core/options.py in render(cls, obj)
   1403         data, metadata = {}, {}
   1404         for hook in hooks:
-> 1405             ret = hook(obj)
   1406             if ret is None:
   1407                 continue

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in pprint_display(obj)
    280     if not ip.display_formatter.formatters['text/plain'].pprint:
    281         return None
--> 282     return display(obj, raw_output=True)
    283 
    284 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in display(obj, raw_output, **kwargs)
    256     elif isinstance(obj, (HoloMap, DynamicMap)):
    257         with option_state(obj):
--> 258             output = map_display(obj)
    259     elif isinstance(obj, Plot):
    260         output = render(obj)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in wrapped(element)
    144         try:
    145             max_frames = OutputSettings.options['max_frames']
--> 146             mimebundle = fn(element, max_frames=max_frames)
    147             if mimebundle is None:
    148                 return {}, {}

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in map_display(vmap, max_frames)
    204         return None
    205 
--> 206     return render(vmap)
    207 
    208 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in render(obj, **kwargs)
     66         renderer = renderer.instance(fig='png')
     67 
---> 68     return renderer.components(obj, **kwargs)
     69 
     70 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/renderer.py in components(self, obj, fmt, comm, **kwargs)
    408                 doc = Document()
    409                 with config.set(embed=embed):
--> 410                     model = plot.layout._render_model(doc, comm)
    411                 if embed:
    412                     return render_model(model, comm)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/viewable.py in _render_model(self, doc, comm)
    453         if comm is None:
    454             comm = state._comm_manager.get_server_comm()
--> 455         model = self.get_root(doc, comm)
    456 
    457         if config.embed:

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/viewable.py in get_root(self, doc, comm, preprocess)
    510         """
    511         doc = init_doc(doc)
--> 512         root = self._get_model(doc, comm=comm)
    513         if preprocess:
    514             self._preprocess(root)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/layout/base.py in _get_model(self, doc, root, parent, comm)
    120         if root is None:
    121             root = model
--> 122         objects = self._get_objects(model, [], doc, root, comm)
    123         props = dict(self._init_params(), objects=objects)
    124         model.update(**self._process_param_change(props))

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/layout/base.py in _get_objects(self, model, old_objects, doc, root, comm)
    110             else:
    111                 try:
--> 112                     child = pane._get_model(doc, root, model, comm)
    113                 except RerenderError:
    114                     return self._get_objects(model, current_objects[:i], doc, root, comm)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/pane/holoviews.py in _get_model(self, doc, root, parent, comm)
    261                   if p in Layoutable.param and p != 'name'}
    262         child_pane = self._panes.get(backend, Pane)(state, **kwargs)
--> 263         self._update_plot(plot, child_pane)
    264         model = child_pane._get_model(doc, root, parent, comm)
    265         if ref in self._plots:

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/pane/holoviews.py in _update_plot(self, plot, pane)
    197             if plot.comm or state._unblocked(plot.document):
    198                 with unlocked():
--> 199                     plot.update(key)
    200                 if plot.comm and 'embedded' not in plot.root.tags:
    201                     plot.push()

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/plot.py in update(self, key)
    981         if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
    982             return self.initialize_plot()
--> 983         item = self.__getitem__(key)
    984         self.traverse(lambda x: setattr(x, '_updated', True))
    985         return item

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/plot.py in __getitem__(self, frame)
    444         if not isinstance(frame, tuple):
    445             frame = self.keys[frame]
--> 446         self.update_frame(frame)
    447         return self.state
    448 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in update_frame(self, key, ranges, element)
   2452                     not subplot in self.dynamic_subplots):
   2453                     continue
-> 2454             subplot.update_frame(key, ranges, element=el)
   2455 
   2456         if not self.batched and isinstance(self.hmap, DynamicMap) and items:

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in update_frame(self, key, ranges, plot, element)
   1541             self._update_hover(element)
   1542 
-> 1543         self._update_glyphs(element, ranges, self.style[self.cyclic_index])
   1544         self._execute_hooks(element)
   1545 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in _update_glyphs(self, element, ranges, style)
   1448             data, mapping, style = self.get_batched_data(element, ranges)
   1449         else:
-> 1450             data, mapping, style = self.get_data(element, ranges, style)
   1451 
   1452         # Include old data if source static

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/geoviews/plotting/bokeh/plot.py in get_data(self, element, ranges, style)
    171         if self._project_operation and self.geographic:
    172             element = self._project_operation(element, projection=self.projection)
--> 173         return super(GeoPlot, self).get_data(element, ranges, style)
    174 
    175 

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/raster.py in get_data(self, element, ranges, style)
    270                         yc.append(ys.mean())
    271                 else:
--> 272                     mask.append(False)
    273 
    274             data = {'xs': XS, 'ys': YS, z.name: zvals[np.array(mask)]}

KeyboardInterrupt: 
print(interactive_winter_mean_maps.__doc__) # Print docstring
 Generate interactive maps of winter mean data 
    Note: this function builds off the functions get_winter_data and interactiveArcticMaps.
    
    Args: 
        da (xr.Dataset or xr.DataArray): data to restrict by time; must contain "time" as a coordinate 
        years (list of str): years over which to compute mean (default to unique years in the dataset)
        start_month (str, optional): first month in winter (default to September)
        end_month (str, optional): second month in winter; this is the following calender year after start_month (default to April)
        force_complete_season (bool, optional): require that winter season returns data if and only if all months have data? i.e. if Sep and Oct have no data, return nothing even if Nov-Apr have data? (default to False) 
        clabel (str, optional): colorbar label (default to "long_name" and "units" if given in attributes of da)
        cmap (str, optional): matplotlib colormap to use (default to "viridis")
        colorbar (bool, optional): show colorbar? (default to True)
        vmin (float, optional): minimum on colorbar (default to 0)
        vmax (float, optional): maximum on colorbar (default to 4)
        title (str, optional): main title to give plot (default to no title)
        ylim (tuple, optional): limits of yaxis in the form min latitude, max latitude (default to (60,90))
        frame_width (int, optional): width of frame. sets figure size of each map (default to 250)
        slider (bool, optional): if da has more than one time coordinate, display maps with a slider? (default to True)
        cols (int, optional): how many columns to show before wrapping, if da has more than one time coordinate (default to 3)
    
    Returns: 
        pl_means (Holoviews map)
    
    
interactive_winter_mean_maps(book_ds_innerArctic["ice_thickness_smoothed"], 
                             years=years,  
                             slider=True, 
                             frame_width=500, 
                             title="ICESat-2 sea ice thickness winter means")

4) Differences in winter sea ice thickness

Here, we compare the sea ice thickness across winter seasons and generate a map of the gridcells differences to identidy any regions of particularly high differences in either direction. This allows us to get an idea of how and where two winter seasons differed in thickness.

Because ICESat-2 didn’t start collecting data until November 2018, we’ll show means for November through April, instead of Sep through April, which we would normally define as the winter season.

The gridcell difference map shows the middle plot minus the first plot, so a positive (colored red in the map) difference would indicate that winter 2019-2020 had a higher sea ice thickness value than winter 2018-2019, and a negative (blue) difference would indicate a lower sea ice thickness value.

is2_winter_means = compute_gridcell_winter_means(book_ds_innerArctic["ice_thickness_smoothed"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds_innerArctic["piomas_ice_thickness"], years=years, start_month="Nov")

for i in range(len(years)): 
    data1 = is2_winter_means.isel(time=i)
    data2 = pio_winter_means.isel(time=i)
    pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
    pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item()) 
    diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference") 
    data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
    data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
    display(data_pl)
    print("\n")



5) Mean monthly sea ice thickness

Compare the monthly means across different winter seasons for the uninterpolated ICESat-2 sea ice thickness data

print(interactive_winter_comparison_lineplot.__doc__)
 Make a bokeh lineplot with markers comparing monthly mean data across winter seasons 
    
    Args: 
        da (xr.DataArray): data to plot and compute mean for; must contain "time" as a coordinate 
        years (list of str): list of years for which to plot data. 2020 would correspond to the winter season defined by start month 2020 - end month 2021 (default to all unique years in da)
        title (str, optional): title to give plot (default to "Winter comparison") 
        frame_width (int, optional): width of figure (default to 600) 
        frame_height (int, optional): height of figure (default to 350) 
        start_month (str, optional): first month in winter (default to September)
        end_month (str, optional): second month in winter; this is the following calender year after start_month (default to April)
        force_complete_season (bool, optional): require that winter season returns data if and only if all months have data? i.e. if Sep and Oct have no data, return nothing even if Nov-Apr have data? (default to False) 
        
       Returns: 
           pl (bokeh lineplot) 
        
    
pl = interactive_winter_comparison_lineplot(book_ds_innerArctic.ice_thickness, 
                            years=[2018,2019,2020],
                            title="ICESat-2 sea ice thickness winter comparison")
display(pl)